require(plm)
require(runjags)
require(coda)
require(wfe)

#### data preparation ####
DPJ.data.Study1 <- read.csv("DPJ_data_Study1.csv")

# exclude periods with no speech
# PS = plenary session, BC = Budget Committee, OC = other committees
DPJ.personnel.period.date.Study1 <- as.Date(c("2000-06-25", "2000-09-09", "2002-09-23", 
                                              "2002-12-10", "2003-12-11", "2004-05-18", 
                                              "2004-09-13", "2005-09-17", "2006-04-07", 
                                              "2006-09-25", "2007-08-31", "2008-09-21", 
                                              "2009-05-16", "2009-09-16", "2010-06-08", 
                                              "2010-09-17", "2011-01-14", "2011-09-02", 
                                              "2012-01-13", "2012-06-04", "2012-10-01", 
                                              "2012-12-26", "2013-09-04", "2014-09-16", 
                                              "2015-01-18", "2015-10-13", "2016-03-27"))

DPJ.PS.data.Study1 <- DPJ.BC.data.Study1 <- DPJ.OC.data.Study1 <- DPJ.data.Study1
for (i in 1:(length(DPJ.personnel.period.date.Study1) - 1)) {
  temp.data <- subset(DPJ.data.Study1, time == i)
  if (sum(temp.data$PS.char > 0) == 0) {
    DPJ.PS.data.Study1 <- subset(DPJ.PS.data.Study1, time != i)
  }
  if (sum(temp.data$BC.char > 0) == 0) {
    DPJ.BC.data.Study1 <- subset(DPJ.BC.data.Study1, time != i)
  }
  if (sum(temp.data$OC.char > 0) == 0) {
    DPJ.OC.data.Study1 <- subset(DPJ.OC.data.Study1, time != i)
  }
}

# convert data frames to the pdata.frame class of the plm package
DPJ.PS.pdata.Study1 <- pdata.frame(DPJ.PS.data.Study1, c("name_jp", "time"))
DPJ.BC.pdata.Study1 <- pdata.frame(DPJ.BC.data.Study1, c("name_jp", "time"))
DPJ.OC.pdata.Study1 <- pdata.frame(DPJ.OC.data.Study1, c("name_jp", "time"))

#### main analyses (Online Appendix E.2) ####
Study1.DPJ.PS <- plm(PS.char ~ leader + chair + 
                       minister + jr.minister + age + 
                       log.term + vote.share * SMD + zombie + 
                       E.43 + E.44 + E.45 + E.46 + E.47, 
                     data = DPJ.PS.pdata.Study1, 
                     model = "within")
round(summary(Study1.DPJ.PS, 
              vcov = vcovHC(Study1.DPJ.PS, 
                            method = "arellano", 
                            type = "HC2"))$coefficients, 3)
nrow(DPJ.PS.pdata.Study1)  # number of observations
length(unique(DPJ.PS.pdata.Study1$name_jp))  # number of MPs
length(unique(DPJ.PS.pdata.Study1$time))  # number of periods

round(exp(Study1.DPJ.PS$coefficients[1]), 3)  # substantive interpretation

#### analysis of committees (Online Appendix E.3) ####
## Budget Committee
Study1.DPJ.BC <- plm(BC.char ~ leader + chair + 
                       minister + jr.minister + age + 
                       log.term + vote.share * SMD + zombie + 
                       E.43 + E.44 + E.45 + E.46 + E.47, 
                     data = DPJ.BC.pdata.Study1, 
                     model = "within")
round(summary(Study1.DPJ.BC, 
              vcov = vcovHC(Study1.DPJ.BC, 
                            method = "arellano", 
                            type = "HC2"))$coefficients, 3)
nrow(DPJ.BC.pdata.Study1)  # number of observations
length(unique(DPJ.BC.pdata.Study1$name_jp))  # number of MPs
length(unique(DPJ.BC.pdata.Study1$time))  # number of periods

round(exp(Study1.DPJ.BC$coefficients[1]), 3)  # substantive interpretation

## other committees
Study1.DPJ.OC <- plm(OC.char ~ leader + chair + 
                       minister + jr.minister + age + 
                       log.term + vote.share * SMD + zombie + 
                       E.43 + E.44 + E.45 + E.46 + E.47, 
                     data = DPJ.OC.pdata.Study1, 
                     model = "within")
round(summary(Study1.DPJ.OC, 
              vcov = vcovHC(Study1.DPJ.OC, 
                            method = "arellano", 
                            type = "HC2"))$coefficients, 3)
nrow(DPJ.OC.pdata.Study1)  # number of observations
length(unique(DPJ.OC.pdata.Study1$name_jp))  # number of MPs
length(unique(DPJ.OC.pdata.Study1$time))  # number of periods

#### robustness check: other dependent variables (Online Appendix E.4) ####
## number of speeches
Study1.DPJ.PS.count <- plm(PS.count ~ leader + chair + 
                             minister + jr.minister + age + 
                             log.term + vote.share * SMD + zombie + 
                             E.43 + E.44 + E.45 + E.46 + E.47, 
                           data = DPJ.PS.pdata.Study1, 
                           model = "within")
round(summary(Study1.DPJ.PS.count, 
              vcov = vcovHC(Study1.DPJ.PS.count, 
                            method = "arellano", 
                            type = "HC2"))$coefficients, 3)

Study1.DPJ.BC.count <- plm(BC.count ~ leader + chair + 
                             minister + jr.minister + age + 
                             log.term + vote.share * SMD + zombie + 
                             E.43 + E.44 + E.45 + E.46 + E.47, 
                           data = DPJ.BC.pdata.Study1, 
                           model = "within")
round(summary(Study1.DPJ.BC.count, 
              vcov = vcovHC(Study1.DPJ.BC.count, 
                            method = "arellano", 
                            type = "HC2"))$coefficients, 3)

Study1.DPJ.OC.count <- plm(OC.count ~ leader + chair + 
                             minister + jr.minister + age + 
                             log.term + vote.share * SMD + zombie + 
                             E.43 + E.44 + E.45 + E.46 + E.47, 
                           data = DPJ.OC.pdata.Study1, 
                           model = "within")
round(summary(Study1.DPJ.OC.count, 
              vcov = vcovHC(Study1.DPJ.OC.count, 
                            method = "arellano", 
                            type = "HC2"))$coefficients, 3)

## making one or more speeches
Study1.DPJ.PS.bin <- plm(PS.bin ~ leader + chair + 
                           minister + jr.minister + age + 
                           log.term + vote.share * SMD + zombie + 
                           E.43 + E.44 + E.45 + E.46 + E.47, 
                         data = DPJ.PS.pdata.Study1, 
                         model = "within")
round(summary(Study1.DPJ.PS.bin, 
              vcov = vcovHC(Study1.DPJ.PS.bin, 
                            method = "arellano", 
                            type = "HC2"))$coefficients, 3)

Study1.DPJ.BC.bin <- plm(BC.bin ~ leader + chair + 
                           minister + jr.minister + age + 
                           log.term + vote.share * SMD + zombie + 
                           E.43 + E.44 + E.45 + E.46 + E.47, 
                         data = DPJ.BC.pdata.Study1, 
                         model = "within")
round(summary(Study1.DPJ.BC.bin, 
              vcov = vcovHC(Study1.DPJ.BC.bin, 
                            method = "arellano", 
                            type = "HC2"))$coefficients, 3)

Study1.DPJ.OC.bin <- plm(OC.bin ~ leader + chair + 
                           minister + jr.minister + age + 
                           log.term + vote.share * SMD + zombie + 
                           E.43 + E.44 + E.45 + E.46 + E.47, 
                         data = DPJ.OC.pdata.Study1, 
                         model = "within")
round(summary(Study1.DPJ.OC.bin, 
              vcov = vcovHC(Study1.DPJ.OC.bin, 
                            method = "arellano", 
                            type = "HC2"))$coefficients, 3)

#### robustness check: censored regression (Online Appendix E.4) ####
# function to draw initial values
initial.fun.Tobit <- function(seed1, seed2, n.var, n.id) {
  set.seed(seed1)
  alpha <- rnorm(1)
  beta <- rnorm(n.var)
  gamma <- rnorm(n.id)
  sigma <- runif(1, 1, 2)
  tau <- runif(1, 1, 2)
  list(alpha = alpha, beta = beta, 
       gamma = gamma, sigma = sigma, tau = tau, 
       .RNG.name = "base::Wichmann-Hill", .RNG.seed = seed2)
}
# draw initial values
Study1.DPJ.PS.Tobit.inits <- list()
for (i in 1:3) {
  Study1.DPJ.PS.Tobit.inits[[i]] <- 
    initial.fun.Tobit(100 * i, 1000 * i, 17, 
                      max(as.numeric(factor(DPJ.PS.data.Study1$name_jp))))
}
# MCMC sampling
Study1.DPJ.PS.Tobit <- run.jags(model = "random_effects_Tobit_JAGS.R", 
                                monitor = c("beta", "alpha", "sigma", "tau"), 
                                data = list(y = ifelse(DPJ.PS.data.Study1$PS.char > 0, DPJ.PS.data.Study1$PS.char, NA), 
                                            z = (DPJ.PS.data.Study1$PS.char > 0) * 1, 
                                            x = scale(as.matrix(DPJ.PS.data.Study1[, c("leader", "chair", "minister", "jr.minister", 
                                                                                       "female", "age", "dynastic", 
                                                                                       "log.term", "vote.share", "SMD", "SMD.VS", "zombie", 
                                                                                       "E.43", "E.44", "E.45", "E.46", "E.47")]), 
                                                      scale = FALSE), 
                                            id = as.numeric(factor(DPJ.PS.data.Study1$name_jp)), 
                                            N = nrow(DPJ.PS.data.Study1), 
                                            n.var = 17, 
                                            n.id = max(as.numeric(factor(DPJ.PS.data.Study1$name_jp)))), 
                                inits = Study1.DPJ.PS.Tobit.inits, 
                                n.chains = 3, burnin = 5000, sample = 5000, modules = "glm", 
                                summarise = FALSE, thin = 1, method = "parallel")
# convergence check
print(which(gelman.diag(Study1.DPJ.PS.Tobit, multivariate = FALSE)[[1]][, 1] > 1.1))
# posterior summary
round(summary(Study1.DPJ.PS.Tobit), 3)
# posterior probability that the coefficient of party leader (beta[1]) is negative
round(mean(unlist(Study1.DPJ.PS.Tobit$mcmc[, 1]) < 0), 3)

Study1.DPJ.BC.Tobit.inits <- list()
for (i in 1:3) {
  Study1.DPJ.BC.Tobit.inits[[i]] <- 
    initial.fun.Tobit(100 * i, 1000 * i, 17, 
                      max(as.numeric(factor(DPJ.BC.data.Study1$name_jp))))
}
Study1.DPJ.BC.Tobit <- run.jags(model = "random_effects_Tobit_JAGS.R", 
                                monitor = c("beta", "alpha", "sigma", "tau"), 
                                data = list(y = ifelse(DPJ.BC.data.Study1$BC.char > 0, DPJ.BC.data.Study1$BC.char, NA), 
                                            z = (DPJ.BC.data.Study1$BC.char > 0) * 1, 
                                            x = scale(as.matrix(DPJ.BC.data.Study1[, c("leader", "chair", "minister", "jr.minister", 
                                                                                       "female", "age", "dynastic", 
                                                                                       "log.term", "vote.share", "SMD", "SMD.VS", "zombie", 
                                                                                       "E.43", "E.44", "E.45", "E.46", "E.47")]), 
                                                      scale = FALSE), 
                                            id = as.numeric(factor(DPJ.BC.data.Study1$name_jp)), 
                                            N = nrow(DPJ.BC.data.Study1), 
                                            n.var = 17, 
                                            n.id = max(as.numeric(factor(DPJ.BC.data.Study1$name_jp)))), 
                                inits = Study1.DPJ.BC.Tobit.inits, 
                                n.chains = 3, burnin = 5000, sample = 5000, modules = "glm", 
                                summarise = FALSE, thin = 1, method = "parallel")
print(which(gelman.diag(Study1.DPJ.BC.Tobit, multivariate = FALSE)[[1]][, 1] > 1.1))
round(summary(Study1.DPJ.BC.Tobit), 3)
round(mean(unlist(Study1.DPJ.BC.Tobit$mcmc[, 1]) < 0), 3)

Study1.DPJ.OC.Tobit.inits <- list()
for (i in 1:3) {
  Study1.DPJ.OC.Tobit.inits[[i]] <- 
    initial.fun.Tobit(100 * i, 1000 * i, 17, 
                      max(as.numeric(factor(DPJ.OC.data.Study1$name_jp))))
}
Study1.DPJ.OC.Tobit <- run.jags(model = "random_effects_Tobit_JAGS.R", 
                                monitor = c("beta", "alpha", "sigma", "tau"), 
                                data = list(y = ifelse(DPJ.OC.data.Study1$OC.char > 0, DPJ.OC.data.Study1$OC.char, NA), 
                                            z = (DPJ.OC.data.Study1$OC.char > 0) * 1, 
                                            x = scale(as.matrix(DPJ.OC.data.Study1[, c("leader", "chair", "minister", "jr.minister", 
                                                                                       "female", "age", "dynastic", 
                                                                                       "log.term", "vote.share", "SMD", "SMD.VS", "zombie", 
                                                                                       "E.43", "E.44", "E.45", "E.46", "E.47")]), 
                                                      scale = FALSE), 
                                            id = as.numeric(factor(DPJ.OC.data.Study1$name_jp)), 
                                            N = nrow(DPJ.OC.data.Study1), 
                                            n.var = 17, 
                                            n.id = max(as.numeric(factor(DPJ.OC.data.Study1$name_jp)))), 
                                inits = Study1.DPJ.OC.Tobit.inits, 
                                n.chains = 3, burnin = 5000, sample = 5000, modules = "glm", 
                                summarise = FALSE, thin = 1, method = "parallel")
print(which(gelman.diag(Study1.DPJ.OC.Tobit, multivariate = FALSE)[[1]][, 1] > 1.1))
round(summary(Study1.DPJ.OC.Tobit), 3)
round(mean(unlist(Study1.DPJ.OC.Tobit$mcmc[, 1]) < 0), 3)

#### robustness check: weighted linear fixed-effect model (Online Appendix E.4) ####
Study1.DPJ.PS.WFE <- wfe(PS.char ~ leader + chair + 
                           minister + jr.minister + age + 
                           log.term + vote.share * SMD + zombie + 
                           E.43 + E.44 + E.45 + E.46 + E.47, 
                         data = DPJ.PS.data.Study1, 
                         unit.index = "name_jp", time.index = "time", 
                         treat = "leader", estimator = "fd")
round(summary(Study1.DPJ.PS.WFE)$coefficients, 3)

Study1.DPJ.BC.WFE <- wfe(BC.char ~ leader + chair + 
                           minister + jr.minister + age + 
                           log.term + vote.share * SMD + zombie + 
                           E.43 + E.44 + E.45 + E.46 + E.47, 
                         data = DPJ.BC.data.Study1, 
                         unit.index = "name_jp", time.index = "time", 
                         treat = "leader", estimator = "fd")
round(summary(Study1.DPJ.BC.WFE)$coefficients, 3)

Study1.DPJ.OC.WFE <- wfe(OC.char ~ leader + chair + 
                           minister + jr.minister + age + 
                           log.term + vote.share * SMD + zombie + 
                           E.43 + E.44 + E.45 + E.46 + E.47, 
                         data = DPJ.OC.data.Study1, 
                         unit.index = "name_jp", time.index = "time", 
                         treat = "leader", estimator = "fd")
round(summary(Study1.DPJ.OC.WFE)$coefficients, 3)